Image Processing with MATLAB

Lecture 35

Author
Affiliation

Minjae Park

Auburn University
MATH 2660 - Spring 2026

Published

April 13, 2026

Overview

Goal

  • Today we will treat an image as a matrix (or a stack of matrices).
  • We will use MATLAB to apply linear algebra ideas to images:
    • linear transformations
    • Fourier transform
    • singular value decomposition (SVD)
  • The emphasis is not on polished code.
  • The emphasis is on seeing how our abstract ideas become concrete.

Class plan

  1. Load an image and inspect it as data
  2. Manipulate color channels
  3. Apply linear transformations to the coordinates
  4. Compute a 2D Fourier transform and filter frequencies
  5. Compress the image using truncated SVD

Logistics

  • If possible, use your own image on your laptop.
  • If you do not have one ready, use a built-in MATLAB image such as:
    • peppers.png
    • cameraman.tif
    • saturn.png
  • A companion demo script is available in scripts/image_processing_demo.m.
  • Working in pairs is fine.
  • If MATLAB desktop is not available, MATLAB Online is completely fine.

Part 1: Images as matrices

Warm-up

  • A grayscale image is a matrix: \[ A \in \mathbb{R}^{m\times n}. \]
  • A color image is usually an \(m\times n\times 3\) array.
  • The three slices correspond to:
    • red
    • green
    • blue
  • Linear algebra is already built into the data structure.

Step 1: Load an image

clear; close all; clc;

% Option 1: use your own image file
I = imread("my_photo.jpg");

% Option 2: use a built-in sample image
% I = imread("peppers.png");
% I = imread("cameraman.tif");

figure;
imshow(I);
title("Original image");

size(I)
class(I)
  • size(I) tells you the dimensions.
  • For RGB images, you should see something like m x n x 3.
  • For grayscale images, you will see just m x n.

Step 2: Convert to doubles and grayscale

I = im2double(I);   % rescales to numbers between 0 and 1

if size(I,3) == 3
    G = rgb2gray(I);
else
    G = I;
end

figure;
subplot(1,2,1);
imshow(I);
title("Color image");

subplot(1,2,2);
imshow(G);
title("Grayscale image");
  • We will use the grayscale image G for Fourier and SVD.
  • This keeps the linear algebra simpler.

Step 3: Look at the RGB channels

if size(I,3) == 3
    R = I(:,:,1);
    Gch = I(:,:,2);
    B = I(:,:,3);

    figure;
    subplot(2,2,1); imshow(I);   title("Original");
    subplot(2,2,2); imshow(R);   title("Red channel");
    subplot(2,2,3); imshow(Gch); title("Green channel");
    subplot(2,2,4); imshow(B);   title("Blue channel");
end

Quick experiment

J = I;
J(:,:,1) = 0;   % remove red

K = I;
K(:,:,3) = 0;   % remove blue

figure;
subplot(1,3,1); imshow(I); title("Original");
subplot(1,3,2); imshow(J); title("No red");
subplot(1,3,3); imshow(K); title("No blue");
  • This is a simple but important point:
    • manipulating entries of the array changes the image directly.

Part 2: Linear transformations

Coordinate transformations

  • An image lives on a grid of pixel coordinates.
  • A linear transformation \[ A = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \] acts on those coordinates.
  • Examples:
    • rotation
    • stretching
    • shear
    • reflection
  • The image values are carried along with the coordinates.

Reading the matrix from its columns

  • The first column of \(A\) is \(A\vec{e}_1\).
  • The second column of \(A\) is \(A\vec{e}_2\).
  • So the columns tell us how the horizontal and vertical directions move.
  • Example: \[ A = \begin{bmatrix} 1.15 & 0.35 \\ 0 & 0.85 \end{bmatrix} \] means: \[ A\vec{e}_1 = \langle 1.15,0 \rangle, \qquad A\vec{e}_2 = \langle 0.35,0.85 \rangle. \]
  • That already tells us: stretch a bit in the \(x\)-direction and add a shear.

Easy built-in examples

theta = 30;  % degrees

I_rot = imrotate(I, theta, "bilinear", "crop");
I_big = imresize(I, 1.3);
I_wide = imresize(I, [size(I,1), round(1.4 * size(I,2))]);

figure;
subplot(2,2,1); imshow(I);      title("Original");
subplot(2,2,2); imshow(I_rot);  title("Rotated");
subplot(2,2,3); imshow(I_big);  title("Scaled up");
subplot(2,2,4); imshow(I_wide); title("Stretched horizontally");
  • These are convenient, but they hide the linear algebra.
  • Next we will specify the matrix ourselves.

Idea before code

  • We build the new image one output pixel at a time.
  • For each output location \[ q = \langle x_{\text{new}}, y_{\text{new}} \rangle, \] we ask:
    • which point \(p\) in the old image landed here?
  • Since the forward rule is \[ q = Ap, \] we solve backwards: \[ p = A^{-1}q. \]
  • Then we copy the color of the old pixel nearest to \(p\).

Custom linear map with a matrix

[m,n,channels] = size(I);
cx = (n + 1)/2;
cy = (m + 1)/2;

A = [1.15 0.35;
     0.00 0.85];   % stretch + shear

Ainv = inv(A);
J = zeros(size(I));

for row_new = 1:m
    for col_new = 1:n
        q = [col_new - cx;
             row_new - cy];

        p = Ainv * q;

        col_old = round(p(1) + cx);
        row_old = round(p(2) + cy);

        if 1 <= row_old && row_old <= m && 1 <= col_old && col_old <= n
            for c = 1:channels
                J(row_new, col_new, c) = I(row_old, col_old, c);
            end
        end
    end
end

figure;
subplot(1,2,1); imshow(I); title("Original");
subplot(1,2,2); imshow(J); title("Applied matrix A");

Why do we use Ainv?

  • We are using inverse mapping.
  • For each pixel in the output image, we ask:
    • where did it come from in the input image?
  • This avoids holes and blank gaps.
  • In the simple code above, we use round(...), so each new pixel copies the nearest old pixel.

What each line is doing

  • cx, cy: move the origin to the center of the image
  • q: one output pixel written as a centered coordinate vector
  • p = Ainv*q: trace that output point back to the source image
  • round(...): choose the nearest source pixel
  • J(row_new, col_new, :) = I(row_old, col_old, :): copy the color

Try one of these matrices

% Rotation by angle theta
theta = pi/6;
A = [cos(theta) -sin(theta);
     sin(theta)  cos(theta)];

% Stretch
% A = [1.4 0.0;
%      0.0 0.7];

% Shear
% A = [1.0 0.4;
%      0.0 1.0];

% Reflection across x-axis
% A = [1  0;
%      0 -1];

Questions

  • Which transformations preserve shape best?
  • Which ones preserve area?
  • Which ones noticeably change orientation?
  • What happens when \(\det(A)=0\) or is very close to \(0\)?

Part 3: Fourier transform

Frequency viewpoint

  • The Fourier transform rewrites the image as a sum of waves.
  • Low frequencies capture broad structure.
  • High frequencies capture rapid changes such as edges, texture, and noise.
  • For a grayscale image matrix \(G\), we use \[ F = \operatorname{fft2}(G). \]

What fft2 is doing

  • The original image is stored in the standard pixel basis.
  • The Fourier transform changes to a different basis:
    • horizontal oscillations
    • vertical oscillations
    • diagonal oscillations
  • So each entry of \(F\) tells us:
    • which frequency is present
    • how strongly it appears
  • In other words, fft2 does for images what Fourier series do for 1D signals.

What we read from the spectrum

  • abs(F) measures the strength of each frequency.
  • angle(F) measures the phase.
  • For our purposes, magnitude is the easiest thing to visualize first.
  • After fftshift(...), the center of the picture corresponds to low frequencies.
  • Frequencies farther away from the center correspond to more rapid oscillation.

Visualizing the spectrum

F = fftshift(fft2(G));
S = log(1 + abs(F));

figure;
subplot(1,2,1);
imshow(G, []);
title("Grayscale image");

subplot(1,2,2);
imshow(S, []);
title("Log magnitude of Fourier spectrum");
colormap parula;
colorbar;
  • The center of the shifted spectrum corresponds to low frequencies.
  • Points farther away correspond to higher frequencies.
  • We use log(1 + abs(F)) because the raw coefficients can vary a lot in size.
  • Without the log scale, a few very large coefficients can hide everything else.

Low-pass vs high-pass

  • A low-pass filter keeps coefficients near the center.
  • So it keeps slowly varying information:
    • overall shape
    • lighting
    • large smooth regions
  • A high-pass filter removes that coarse information and keeps rapid changes:
    • edges
    • texture
    • noise
  • This is why low-pass filtering blurs, while high-pass filtering highlights outlines.

Low-pass filter

[m,n] = size(G);
[u,v] = meshgrid(-floor(n/2):ceil(n/2)-1, -floor(m/2):ceil(m/2)-1);
R = sqrt(u.^2 + v.^2);

D0 = 40;              % cutoff radius
Hlow = double(R <= D0);

Flow = F .* Hlow;
Glow = real(ifft2(ifftshift(Flow)));

figure;
subplot(1,3,1); imshow(G, []);     title("Original");
subplot(1,3,2); imshow(Hlow, []);  title("Low-pass mask");
subplot(1,3,3); imshow(Glow, []);  title("Low-pass result");

High-pass filter

Hhigh = double(R >= D0);
Fhigh = F .* Hhigh;
Ghigh = real(ifft2(ifftshift(Fhigh)));

figure;
subplot(1,3,1); imshow(G, []);      title("Original");
subplot(1,3,2); imshow(Hhigh, []);  title("High-pass mask");
subplot(1,3,3); imshow(mat2gray(Ghigh)); title("High-pass result");
  • Low-pass filtering tends to blur the image.
  • High-pass filtering emphasizes edges and fine details.

Questions

  • What visual information is lost under low-pass filtering?
  • What kinds of features survive high-pass filtering?
  • How does this compare to keeping only a few large singular values?

Part 4: SVD compression

Truncated SVD

  • For a grayscale image matrix \(G\), compute \[ G = U\Sigma V^T. \]
  • Keeping only the first \(r\) singular values gives \[ G_r = U_r \Sigma_r V_r^T. \]
  • This is the best rank-\(r\) approximation in least-squares sense.

How to read the SVD for an image

  • Write the SVD as \[ G = \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T + \cdots + \sigma_r u_r v_r^T. \]
  • Each term \(\sigma_i u_i v_i^T\) is a rank-one image.
  • You can think of it as:
    • one vertical pattern \(u_i\)
    • one horizontal pattern \(v_i\)
    • multiplied together
  • The singular value \(\sigma_i\) tells us how important that pattern is.

Why the first few singular values matter

  • The singular values are ordered: \[ \sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_r \ge 0. \]
  • So the first few terms capture the dominant large-scale structure.
  • Later terms often carry finer detail, small texture, or noise.
  • Truncating the SVD means:
    • keep the most important patterns
    • throw away weaker ones
  • That gives compression.

Fourier and SVD are different viewpoints

  • Fourier uses a fixed basis of waves.
  • SVD uses a data-adaptive basis built from the image itself.
  • Fourier asks:
    • which frequencies are present?
  • SVD asks:
    • which rank-one patterns best explain the image?
  • Both compress information, but they keep different kinds of structure.

Basic MATLAB code

[U,S,V] = svd(G, "econ");
s = diag(S);

figure;
plot(s / s(1), "LineWidth", 2);
title("Normalized singular values");
xlabel("Index");
ylabel("\sigma_i / \sigma_1");
grid on;

Reconstruct with several ranks

ranks = [5 20 50 100];

figure;
subplot(2,3,1);
imshow(G, []);
title("Original");

for j = 1:length(ranks)
    r = ranks(j);
    Gr = U(:,1:r) * S(1:r,1:r) * V(:,1:r)';
    subplot(2,3,j+1);
    imshow(Gr, []);
    title(sprintf("rank = %d", r));
end

Energy captured

energy = cumsum(s.^2) / sum(s.^2);

figure;
plot(energy, "LineWidth", 2);
xlabel("rank r");
ylabel("cumulative energy");
title("How much of the image is captured?");
grid on;
r = 50;
fprintf("Energy captured by rank-%d approximation: %.4f\n", ...
        r, energy(r));

Why this is called the “best” rank-\(r\) approximation

  • Among all matrices of rank at most \(r\), the truncated SVD approximation \[ G_r = U_r \Sigma_r V_r^T \] is the closest to \(G\) in least-squares sense.
  • So if we insist on using only \(r\) rank-one pieces, SVD is optimal.
  • That is why it is so important in compression and dimension reduction.

Optional: color SVD

R = I(:,:,1); Gch = I(:,:,2); B = I(:,:,3);
r = 40;

[Ur,Sr,Vr] = svd(R, "econ");
[Ug,Sg,Vg] = svd(Gch, "econ");
[Ub,Sb,Vb] = svd(B, "econ");

Rr = Ur(:,1:r) * Sr(1:r,1:r) * Vr(:,1:r)';
Gr = Ug(:,1:r) * Sg(1:r,1:r) * Vg(:,1:r)';
Br = Ub(:,1:r) * Sb(1:r,1:r) * Vb(:,1:r)';

I_svd = cat(3, Rr, Gr, Br);

figure;
subplot(1,2,1); imshow(I);      title("Original");
subplot(1,2,2); imshow(I_svd);  title("Color SVD approximation");

In-class tasks

Core tasks

  1. Load one image and display its size
  2. Convert it to grayscale
  3. Remove one RGB channel or swap two channels
  4. Apply one custom linear transformation matrix
  5. Apply one low-pass or high-pass Fourier filter
  6. Build one rank-\(r\) SVD approximation

Further questions

  • Compare two different images: which one needs a larger rank for similar quality?
  • Try several cutoff radii in the Fourier filter.
  • Compare a low-pass Fourier approximation with a low-rank SVD approximation.
  • Try to explain visually:
    • what Fourier keeps well
    • what SVD keeps well

Some answers

  • Fourier vs SVD?
    • Fourier decomposes into frequencies.
    • SVD decomposes into rank-one patterns adapted to the image.
  • Why does low-pass filtering blur the image?
    • Because sharp edges need high frequencies, and low-pass removes them.
  • Why does high-pass filtering highlight edges?
    • Because edges are places where pixel values change rapidly.
  • Why does truncated SVD still look recognizable?
    • Because the first few singular values often capture most of the important structure.
  • Is Fourier better than SVD?
    • Not universally. They are useful for different goals and preserve different features.
  • Why did we switch to grayscale first?
    • It keeps the mathematics cleaner. For color, we can repeat the same idea channel by channel.

Wrap-up

  • An image can be understood as:
    • a matrix of intensities
    • a function on a grid
    • a sum of frequency modes
    • a sum of rank-one patterns
  • This is why image processing is such a natural playground for linear algebra.